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We present a formulation of the Hartree-Fock-Bogoliubov (HFB) equations 
which solves the problem directly in the basis of natural orbitals. This pro- 
vides a very efficient scheme which is particularly suited for large scale calcu- 
lations on coordinate-space grids. 
The production of new nuclei towards the drip lines and in the super-heavy region (for a 
review see Jl|) has raised a growing interest in the refinement of nuclear mean-field models 
to accommodate the larger body of experimental data. At the level of precision reached 
nowadays, a correct treatment of pairing becomes a crucial ingredient 0. This calls for a 
full variational optimization of the pairing wave function, what is known as the Hartree-Fock- 
Bogoliubov (HFB) treatment 0,|J. It is the aim of this short note to present an efficient 
solution scheme for the HFB equations which relies on a direct variational optimization 
within the natural orbital basis. The following discussion is formulated for one sort of 
Fermions. Nuclear applications will use that scheme for neutrons and protons separately. A 
generalization to proton-neutron pairing is obvious. 

Starting point is the BCS ansatz for the wave function of a pairing many-body system 
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where a+ generates a particle in the state ip a , a is the time reversed partner of a, v a the 
occupation amplitude for the state, and u a = <Jl — v a the non-occupation amplitude. The 
ansatz ([I]) requires that the single-particle states are orthonormalized, 

^ > ( 2 ) 

and that the occupations add up to the total particle number J2 v a = N. These two 
presuppositions will have to be added as boundary conditions in the variation later on. 
Note, furthermore, that we have assumed here the case of a time-even state (ground state 
of even-even nuclei) for which u a and v a can be chosen real and for which we can construct 
the time-reversed wavefunction as tp a = —ia y (p* a where a y is a Pauli spin matrix. 

The BCS ansatz (JJ) carries only one-body information which is summarized in the two 
key densities: the one-body density operator 

P = Y, V l ( P*Vt , ( 3 ) 

and the pair-density operator x — J2 a u a v a L Pa L Pa fr° m which we are going to use only the 
local part 

Xfr) = J2 U a V afi( r )fa( r ) ■ ( 4 ) 

Thereby, we have employed time- reversal symmetry which also yields the property x*( r ) = 
x(r). A discussion of the physical content of the pair density can be found in ||. Several 
useful properties are further discussed in Note that the one-body density matrix @ 
is diagonal in the chosen representation which means that we are dealing with the natural 
orbital basis. It is noteworthy that the pair density is also diagonal in the same basis. This 
is guaranteed by the relation [p, x] = which can be derived on general grounds [|J . 

The standard solution scheme for the HFB problem deals with a four times as large super- 
density-matrix which encompasses p as well as x, and it proceeds in three subsequent unitary 
transformations 0. There is a treatment in terms of wave functions which is particularly 
suited for coordinate space representations and which uses a double set of single-particle 
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wavef unctions, one for the occupation part and one for the non-occupation part Both 
schemes are plagued by the doubling of the representation which adds substantial overload to 
the calculations and which becomes very cumbersome in deformed nuclei. The BCS-ansatz 
(|I|) can be formulated in terms of a single set of wave functions, the natural orbitals, and a 
few occupation amplitudes. It is possible to keep the scheme at that level of expense by a 
direct variational exploitation of the ansatz, as will be shown in the following. 

For simplicity, we work here with the case that the energy separates into a mean-field 
part and a pairing part as 

E = E m{ (p) + E paiT {x) . (5) 

Although the present scheme is applicable under more general conditions, we use now a 
particular model for the pairing energy, a volume pairing with a zero-range force 

£pair = V P I d 3 r X (r) 2 (6) 
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which employs only the local pair density (|j). The energy functional determines the mean- 
field Hamiltonian and pair potential as the functional derivatives 

f dE .. , dE 1 . . 

""'=W ' A(r) = B^) = 2 VpXir) ' (?) 
Note that we obtain a local pair potential for the particular pairing functional @. 

The variation with respect to the occupation amplitudes v a yields 
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where ep is the Lagrange multiplier for the particle number constraint. This equation can 
be resolved in the standard manner and yields 
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- T - , Ka ~ €¥ . (8) 



Note that the gap potential A does not necessarily commute with the mean-field Hamiltonian 
h. Only the diagonal elements in the natural orbital basis enter and no information about 
the non-diagonal elements is needed. 



The variation with respect to the single particle wavefunctions needs to take care of 
the orthonormality (|2|) which is done by adding - E a /3 A aj3 {^P a \fp)- The thus constrained 
variation yields 

ti a ¥a = E/J^a/jV/J (9) 

with a generalized mean-field Hamiltonian 

dE r - i 
7-L a ¥> a = g—^: = [v 2 a h m{ + u a v a A(r)\ip a . (10) 

This TC a is a state dependent Hamiltonian and the full matrix of Lagrangian multipliers 
needs to be taken into account. Thereby it is crucial that they constitute a symmetric 
matrix X a/3 = \a a . This allows to symmetrize explicitely 



U^\n a + n^ a ) (ii) 



Kp 2 

which links pairwise all a with f3 and thus eq. (0) delivers a decisive problem. 

Altogether, the equations (||), fliTf), ([]), and (jTT|) constitute the HFB equations formu- 
lated directly in the natural orbital basis. This nonlinear problem is best solved iteratively. 
We propose an interlaced iteration employing the efficient damped gradient step which is 
best suited for coordinate space techniques 0: 

1. Start from a given set {<£ a v a }. 

2. Compute the densities p and x( r )> an d subsequently the corresponding Hamiltonians 
h m { and A(r). 

3. Compute the new v a and u a according to to Eq. (^). 

4. Compute the action of the state-dependent mean field 1rt a and store the resulting set 

5. Compute the matrix of constraints X a/3 , Eq. (|TTJ) . 
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6. Perform the damped gradient step with 

v a - ow a - v a \n aVa - E/3 KpVp]} (12) 

where O means orthonormalization and T> a is an appropriate damping operator. This 
completes the scheme and returns to the starting point, step ffl. 



The state-dependent Hamiltonian (|10[) requires a state-dependent damping for which we 
generalize the form of || to 

X> = ^ (13) 

° v 2 a (50 MeV + f) + M a t;4max{A(r)} 

where T is the kinetic energy operator and xq ~0.2 a numerical parameter. With that choice, 
we have implemented this scheme successfully into a spherical Skyrme-Hartree-Fock code 
and tested it extensively for a variety of Skyrme forces and nuclei from 16 O to the isotopes 
of Pb, including proton rich as well as neutron rich exotic nuclei. The scheme proves to be 
reliable and very efficient. It allows a fast computation of the HFB ground state, costing 
only 20% more iterations than the much simpler BCS approximation, and each iteration as 
such is as fast as in the BCS case because only one set of single particle wavefunctions is 
handled. We thus have a promising alternative HFB scheme which can simplify large scale 
calculations of deformed nuclei. 
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